Metamagnetism in the 2D Hubbard Model with easy axis 
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Although the Hubbard model is widely investigated, there are surprisingly few attempts to study 
the behavior of such a model in an external magnetic field. Using the Projector Quantum Monte 
Carlo technique, we show that the Hubbard model with an easy axis exhibits metamagnetic behavior 
if an external field is turned on. For the case of intermediate correlations strength U, we observe a 
smooth transition from an antiferromagnetic regime to a paramagnetic phase. While the staggered 
magnetization will decrease linearly up to a critical field Be, uniform magnetization develops only 
for fields higher than Be- 
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I. INTRODUCTION 



The Hubbard model is one-of the simplest models of 
strongly correlated electrons.!!! The magnetic properties 
of this. |model have been extensively studied for many 
years.!3Q But only in a few instances the influence of an 
external magnetic field Jaeing coupled to the electrons, 
has been investigated.ElQ A very popular approach is the 
Peierls substitution, i.e. a hopping amplitude of the elec- 
trons which depends on the vector potential of the exter- 
nal field. This is used, e.g., to study the superconducting 
properties-pf a Hubbard ring or torus threaded by a mag- 
netic Aux.e! It would also be appropriate to calculate Hall 
coefficients in such systems. 

A different approach is to include a Zeeman term in the 
Hamiltonian, i.e. to couple the extesaal magnetic field 
directly to the spins of the electrons .!a This case is well 
suited for calculating static properties, such as magneti- 
zation. 

Since many years it is well known that in alloys with a 
layered structure the magnetization shows a specific be- 
havior. If the planes are itself ferromagnetically ordered 
but the coupling between them is antiferromagnetic, one 
observes that in an external field the total magnetization 
first slowly increases linearly, tben suddenly strongly rises 
before saturation takes place.H This was first observed 
by Becquerel andi-jvan den Handel who coined the term 
'met amagnet ism ' 

Especially, since metamagnetic behavior was found in 
heavy fermion compounds, the term 'metamagnetism' is 
used whenever the magnetic susceptibility Xm{B) has a 
maximum at a critical field Be, i.e. the magnetization 
M{B) has a point of infiexion at that field value, even if 
no phase transition occurs. 

It is widely believed that antiferromagnetic correla- 
tions play a crucial role in metamagnetic behavior. Since 
antiferromagnetic correlations are inherent in the Hub- 
bard model, it is a very interesting question to study 
whether the Hubbard model shows metamagnetic behav- 



ior or not. ii 

Recently Held et alxi investigated an anisotropic Hub- 
bard model in a magnetic field coupled via a Zeeman 
interaction term. Using the Grand Canonical Quantum 
Monte Carlo approach, they calculated in d = cxd a mag- 
netic phase diagram and found several phase transitions 
of first and second order. It is an open question whether 
these phases still exist in the more realistic case oi d — 2 
or 3. 

Here, we consider the two-dimensional Hubbard model 
on a square lattice in an external magnetic field B cou- 
pled to the spins of the electrons via a Zeeman term. The 
Hamiltonian H is given by 

U ■ 



H — tijcl^Cj^ + Tlia- 



(1) 



where tij denotes nearest neighbor hopping, B^ is the 
magnetic field parallel to the z axis, and si = <Jni„ 
is the spin in z direction. While the Hamiltonian itself is 
isotropic, an easy axis along the z direction will be intro- 
duced by the simulational procedure as will be discussed 
later on. 



II. METHOD 

Here we briefly review the Projector Quantum Monte 
Carlo (PQMC) method for fermions in the ground-state. 
For a detailed discussion, the reader is referred to Ref. [lC| . 

The key idea of the PQMC algorithm is to project out 
the ground-state wave function | ^E'o) of a lattice fermion 
Hamiltonian H from a given trial wave function j $t) by 
applying the operator exp{—f3H) on | ^t) according to 
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The expectation values of physical quantities A are then 
obtained from 



(A) = lim 

p— >oo 



($^|e-'3«Ae-'3«|$T) 



(3) 



($T|e-2/3J^ 1$^} 

Applying the Trotter-Suzuki decompositionQ^E — and 
the discrete Hubbard-Stratonovich transforniationll3 to 
the projection operator, the effect of the projection op- 
erator on the trial state can be rewritten symbolically 
as a sum over the Hubbard-Stratonovich spins {u}, 
Q-PH I $t) = E{<t} ^iW}) I The expectation value 
of a physical quantity A is then obtained from 

E (a>T|F({a})AF({a'})|$T) 

(A) = i'^i^i'^'y (A] 

^ ' E ($T|F(W)F({a'})|$T) ■ ^ ' 

{'^},{'^'} 

TO|-evaluate these sums, the Monte-Carlo method is 
used,Ej utilizing 

^iW}, W'}) I = I {'^t\ Fi{a}) Fi{a'}) \ $t) | (5) 

as the weight of a configuration of Hubbard-Stratonovich 
spins. Since in general u;{{a}, {cr'}) can be negative for 
some spin configurations {cr}, it can be difficult to eval- 
uate Eq. (jj) numerically. This problem is often referred 
to as the minus-sign problem. 

All Quantum Monte Carlo simulations suffer from the 
so-called 'minus-sign' problem though it does not always 
occur at half-filling. In the PQMC scheme, the minus- 
sign problem can be avoided for the bare Hubbard model 
at half-filling if one uses a spin density wave (SDW) 
ground-state as the trial wave function. In our sim- 
ulations, we found that an appropriately chosen SDW 



ground-state wave function reduces the minus-sign prob- 
lem in case of an additional external magnetic field, too. 

The (zero-field) Hubbard model is invariant under 
SU(2) spin rotation symmetry. Since the zero-field 
Hamiltonian commutes with 5*^, the eigenstates of Sz 
are a natural choice for a basis of states. At half fill- 
ing, the ground-state is a state with 5*2 = 0. Therefore, 
one usually constructs the trial wave function as a direct 
product of spin up and spin down wave functions with an 
equal and fixed number of electrons in each spin direc- 
tion. Hence, the PQMC scheme does not only conserve 
the total number of electrons, but moreover it restricts 
the simulation to states of = 0. 

In order to incorporate the external magnetic field, 
we have to remove this constraint. Therefore, we have 
extended the PQMC algorithm to all eigenstates of S^. 
This is achieved by allowing the number of electrons with 
a certain spin to change while still keeping the total num- 
ber of electrons fixed. Hence, we still work in the canon- 
ical ensemble appropriate for the ground-state. In the 
framework of the PQMC method, our procedure corre- 
sponds to a manifold of trial wave functions, all differing 
in spin Sz, which are all sampled by the Monte Carlo 
method. 

To be more specific, let us write a general trial wave 
function as a sum over trial wave functions with fixed Sz , 



(6) 



Now, the expectation value of an operator A which con- 
serves the spin, reads 
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where the absolute value of 



co{{a}, {a'}, Sz) = [aiSz)f{<i>T{Sz)\F{{a}) Fi{a'}) \ $^(5,)) 



(8) 



is now being used as the generalized weight of a config- 
uration. Application of the Monte-Carlo method is now 
straightforward. We want to stress that each point of 
the configuration space is still characterized by a definite 
value of Sz ■ The original scheme of S'z =0 for a half- filled 
band corresponds to a choice of 

/r, N _ J 1 for 5^ = 
"I zj — 'j^ Q otherwise . ^ ' 



Since the Zeeman term is bilinear in the electronic op- 
erators and commutes with the other parts of the Hamil- 
tonian (^, it is easily incorporated into the operator 
exp(— and the Hubbard-Stratonovich transforma- 
tion stays unchanged. 

The trial wave functions in our scheme are composed 
of direct products of Slater determinants of electrons of 
fixed spin directions, 

|$t(5.)) = |4(5,))®|4(5.)). (10) 
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Hence no linear combinations of up and of down spins 
can occur as would be necessary to construct eigenstates 
of Sx or Sy. This introduces an easy axis along the z axis 
into the simulation, constraining the spins to lie parallel 
to it. Since the Hamiltonian conserves spin directions, 
the structure of the trial wave function also applies to 
the projected ground state wave function, Eq. (^. Con- 
sequently, the easy axis is conserved throughout the sim- 
ulation. 



III. NUMERICAL RESULTS AND DISCUSSION 

We have performed simulations for square lattices up 
to a linear system size of L = 8, i.e. N = L x L lattice 
sites. On average, we used m = 64-128 time slices for 
our Trotter-Suzuki decomposition. The number of elec- 
trons was set to N and kept fixed throughout the whole 
simulation. As outlined above, however, the number of 
electrons with a given spin direction may change during 
the course of the simulation so that a net magnetization 
results. Typically, we run for an initial warm-up and 
following measurements several thousands Monte Carlo 
sweeps. This procedure was repeated about 10 times to 
get independent data from which the average and the 
error was computed. 

In order to compare our results to the work of Held 
et at, we use U — 2 [m units of t). We made extensive 
studies using different projection parameters /3 to ensure 
proper convergence of the energy and the magnetization, 
resp., to the ground-state. It can be observed that the 
energy converges to a /3-independent value much faster 
than the magnetization. However, it turned out that in 
most cases a value of /3 = 6 is sufficient to reach a fi- 
nal value. Furthermore, we checked that the error due 
to the Trotter-Suzuki decomposition is smaller than the 
statistical error in our data. This was achieved by vary- 
ing the number m of time slices. Then the error due to 
the decomposition can be estimated from a scaling of the 
ground-state energy versus Ijw?. 

For the trial state | $t), we choose a spin-density wave 
state. The external magnetic field is parallel to the quan- 
tization axis of the spins of the electrons and thus parallel 
to the easy axis. Applying such an external field to the 
system, the electrons can gain energy by orienting their 
spin in the direction of the field. Therefore, one would ex- 
pect a decrease in the ground-state energy with increas- 
ing field. In Fig. Q we have depicted the ground-state 
energy per site. Eg, as measured for three different sys- 
tem sizes. With increasing magnetic field the energy 
is considerably reduced. For small systems, we observe 
a linear decrease of Eg while for larger systems the slope 
slowly changes. 

Clearly, there are finite size effects for the energy. We 
observe that our data Eg{B, N) scale according to a 1/A^ 
behavior. Extrapolating them to TV — > oo, we derived 
Eg{B) for an infinitely large system, too. For fields of 



B = 0.2t...0.8t, the ground-state energy Eg{B) can 
nicely be fitted with a quadratic function, i.e. Eg{B) cx 
B'^ . This quadratic behavior can easily be understood: 
For free electrons {U = 0), it is known that a Zeeman 
coupling of the electrons leads to a (temperature inde- 
pendent) van Vleck contribution to the static suscepti- 
bility. Therefore, the energy should show a quadratic 
dependence on the magnetic field as confirmed in our 
simulations. This behavior, observed for U = 2t, could 
be a sign that the antiferromagnetic correlations intro- 
duced by the Coulomb repulsion U, are broken up by 
the external field B. So we expect some influence on the 
magnetization as will be discussed below. 
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FIG. 1. Ground-state energy versus external magnetic field 
Bz- The symbols correspond to a linear system size of 
L = 4(x), 6(D), 8(o), and c>o(*). Lines are guides to the 
eye only. 

In our simulation we computed the spin-spin correla- 
tion function 

^(q) = ^E«''''^'''"''^^«"^T -«a)K-T (11) 

In order to extrapola|te to the thermodynamic limit, we 
plot S'(q)/iV vs. l/iV.ml3 It should follow a straight fine 
according to 

^(q) - NTn\ + SM) , (12) 

where Sc is the connected structure factor and rtiq the 
magnetization 

"^« = ^E^"'''''«"*T-"a))- (13) 

i 

From the extrapolated value iV ^ oo, we obtain the 
square of the magnetization m^. We have followed this 
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procedure for g = and q = Q = (tt, tt) to obtain the 
uniform and the staggered magnetization. 

Our results for the magnetizations mo{B) and mQ{B) 
are shown in Fig. ^. 
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FIG. 2. Uniform (x) resp. staggered (□) magnetization 
versus external magnetic field. Note the transition point from 
an antiferromagnetic to a paramagnetic phase at Be ~ 0.25t. 
Lines are guides to the eye only. 

In zero field, the Hubbard model shows antiferromag- 
netic order. With increasing external field, the stag- 
gered magnetization clearly decreases up to a critical field 
Be ~ 0.25i where it vanishes. At about the same field 
value, the uniform magnetization strongly rises. The in- 
flexion point is clearly seen in the uniform magnetization, 
thus a metamagnetic behavior takes place. 

Since in our case the trial wave function is a general- 
ized spin density wave from which the true ground-state 
is projected out, any non-zero staggered magnetization 
means that the system is still in an antiferromagnetic 
ground-state. This is true up to the critical field i?c- 

According to a theorem due to Mermin and Wagner, an 
isotropic two-dimensional system does not undergo a con- 
tinuous transition at any finite temperature. However, 
long-range order at zero temperature is not excluded. Al- 
though the simulation method cannot deal with T = 
directly, the true ground-state is approached for suffi- 
ciently high projection parameters f3. Even if we would 
not have reached high enough (3, the system would behave 
effectively as a long-range ordered one if the correlation 
length is larger than the system size. 

The fact that the system has an easy axis due to the 
simulational constraint, is certainly a limitation. In an 
anisotropic Heisenberg model, being the limiting case of 
large U, there exists a minimum fieLd _Btr at which a 
so-called spin-flop transition occursJia The spins of the 
electrons will then orient themselves perpendicular to the 



external field. In the isotropic model, this will happen 
already for an infinitesimal psjiall field Bt^. Raising the 
anisotropy, i?tr will increase.li3 Further attempts have to 
be made to clarify whether this scenario holds for the 
Hubbard model with intermediate U , too. _ 
If one compares our results with those of Held et al.,U 
there arc remarkable differences. For low temperatures, 
they found at low external fields a constant, finite stag- 
gered magnetization and vanishing uniform magnetiza- 
tion. At a critical field of Be ~ 0.12i, a first order phase 
transition takes place leading to a jump in both magne- 
tization curves. For fields larger than Be, the staggered 
magnetization remains zero while the uniform magnetiza- 
tion increases further. In contrast, our T = value of Be 
is twice as large-as Be and close to the mean-field value 
of w 0.27i.ll3 Besides, we find a rather smooth tran- 
sition in the staggered magnetization decreasing steadily 
up to Be- Due to strong fluctuations and large statistical 
errors close to the phase transition, we were not able to 
resolve the question if there is a mixed phase with mo ^ 
and mg ^ 0. Nevertheless, in our opinion the absence 
of a jump in mg is a strong indication of a second order 
phase transition in two dimensions. 

IV. SUMMARY 

To summarize, we have studied the half-filled two- 
dimensional Hubbard model with an easy axis in an ex- 
ternal magnetic field which was coupled to the spin of the 
electrons via a Zeeman term. The model was investigated 
numerically using an enhanced version of the Projector 
Quantum Monte Carlo method. 

The model shows in zero field an antiferromagnetic 
ground-state which remains present in increasing exter- 
nal magnetic fields up to a critical field value Be « 0.25i. 
For higher fields the system is found to be in a para- 
magnetic state with field-induced spin orientation. Our 
data suggest that the phase transition at Be should be of 
second order. 
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